
# Import logging library and log file prior to initialising rest of code. Typically assists in prevent log file export issues.
# --------------------------------#
import logging
import datetime
logging.basicConfig(filename=datetime.datetime.now().strftime("%Y%m%d_%H%M%S")+'_log_service_hex_joins.log', filemode='w', format='%(asctime)s - %(name)s - %(levelname)s - %(message)s', level=logging.INFO)

# Import remaining libraries.
# --------------------------------#
# General libraries.
import pandas as pd
import numpy as np
from numpy import arange
import json
import time
import sys
import plotly.express as px #fig.show(renderer='notebook') #Allows all plotly plots to be exported with the html.
import missingno as msno
import matplotlib.pyplot as plt
import math

# Spatial libraries.
import geopandas as gpd
import shapely
from shapely.geometry import Polygon, Point, MultiPolygon, MultiPoint, GeometryCollection
from shapely.ops import nearest_points
from shapely.geometry.base import geom_factory
from shapely.geos import lgeos
from shapely import wkt
from area import area

# Mapping.
import folium
from folium import plugins
from folium import Map, FeatureGroup, Marker, LayerControl
from folium.features import GeoJson, GeoJsonTooltip, GeoJsonPopup
from folium.plugins import MarkerCluster
from folium.plugins import MeasureControl
from folium.plugins import MousePosition
from folium.plugins import FloatImage

# Introspection specific libraries.
# --------------------------------#
#Stationary/non-stationary time series analysis.
from statsmodels.tsa.stattools import adfuller

# Ridge Regression.
from sklearn.linear_model import RidgeCV
from sklearn.model_selection import RepeatedKFold

# Auto ARIMA
from statsmodels.tsa.arima_model import ARIMA
import pmdarima as pm #pip install pmdarima --user (resolved environment issues.)

# Classification specific libraries.
# --------------------------------#
from sklearn.cluster import DBSCAN
from sklearn.cluster import KMeans


area_me_rdme={
    'type':float,
    'units':'square kms',
    'docstring':'From geometric object (shapely geometry), in decimal degrees format (crs="EPSG:4326"), outputs area of polygon in square kilometers'
}

def area_me(polygon)->area_me_rdme: #Polygon is in DD format for coords.
    return area(json.dumps(shapely.geometry.mapping(polygon)))/10**6


# To view annotations please uncomment the below eg of docstring below.
# -------------------#

# Export unit format.
# area_me.__annotations__['return']['type']

# Export units.
# area_me.__annotations__['return']['units']

# Export Description.
area_me.__annotations__['return']['docstring']

'From geometric object (shapely geometry), in decimal degrees format (crs="EPSG:4326"), outputs area of polygon in square kilometers'


polyFix_it_rdme={
    'type':'range of folders related to geometric test and manipulation',
    'units':'Varied per column',
    'docstring':'Function which takes in a geometry and geometry index and checks for various geometric errors, applies simple fixes and exports information about the original and adjusted geometries'
}

def polyFix_it(geom, index)->polyFix_it_rdme:
    """
    OUTPUT COLUMNS:
    -------------------------   
    'Orig_Area', - original geometry area in sqkm
    'Orignal Geometry Type', - original geometry type, eg. Point, Polygon etc.
    'Orignal Geometry Validity', - Output of validity (valid/not valid) for original geometry.
    'Orignal Geometry Validity Reason', Validity check eg. self intersection, polygon not fully joined etc.
    'Final Geometry',  - corrected geometry
    'Final Area', - corrected geometry area in sqkm
    'Final Geometry Validity', - Output of validity (valid/not valid) for corrected geometry.
    'Final Valid Ratio' - Ratio of area, original geometry area/corrected geometry area. This should be close to 1.
    
    EXAMPLE OF CALLING FUNCTION:
    ------------------------- 
    # pdf[['Orig_Area', 'Orignal Geometry Type', 'Orignal Geometry Validity', 'Orignal Geometry Validity Reason', 'Final Geometry', 'Final Area', 'Final Geometry Validity', 'Final Valid Ratio']]=pdf.apply(lambda x: polyFix_it(x['geometry'], x.name), axis=1, result_type='expand')
    
    """
    # print(index) #Debug.
    
    # Original area:
    orig_area=area_me(geom)
    orig_area_type=geom.geom_type
    is_valid_=geom.is_valid
    is_valid_reason=lgeos.GEOSisValidReason(geom._geom)
    
    if is_valid_==True:
        new_geom=geom
        corr_method='Valid geometry, no adjustments required.'
        return [orig_area, orig_area_type, is_valid_, is_valid_reason, new_geom, orig_area, 'Valid geometry, no adjustments required.', 1]
    
    # Buffer:
    buff_geom=geom.buffer(0)
    buffer_area=area_me(buff_geom)
    buffer_valid=buff_geom.is_valid
    if buffer_valid==True and buffer_area>0:
        buffer_valid_ratio=round(orig_area/buffer_area,1)
    else:
        buffer_valid_ratio=0
    
    # Stitching:
    if geom.geom_type=='Polygon':
        ls = LineString(geom.exterior)
        lr = LineString(ls.coords[:] + ls.coords[0:1])
        mls = shapely.ops.unary_union(lr)
    #     mls.geom_type
        stitch_geom = MultiPolygon(list(shapely.ops.polygonize(mls)))
        stitch_area=area_me(stitch_geom)
        stitch_valid=stitch_geom.is_valid
        if stitch_valid==True and stitch_area>0:
            stitch_valid_ratio=round(orig_area/stitch_area,1)
        else:
            stitch_valid_ratio=0
        
    #Non-viable option:
    non_viable=[orig_area, orig_area_type, is_valid_, is_valid_reason, 'Geometry Cannot Be Resolved', 0, 'Geometry Cannot Be Resolved', 'Geometry Cannot Be Resolved']
        
    if buffer_area==0 and stitch_area==0:
        return non_viable
    
#     if buffer_area==0: #Debug.
#         print(index)
        
    if orig_area>0 and buffer_area>0:
        buffer_valid_ratio=round(orig_area/buffer_area,1)
        if (is_valid_ == False) and (buffer_valid_ratio<2):
            new_geom=buff_geom
#             corr_method='Buffer used.'
        elif (is_valid_ == False) and (buffer_valid_ratio>=2):
            if geom.geom_type=='Polygon':
                new_geom=stitch_geom
                corr_method='"Stiching Points."'
            else:
                return non_viable
    elif orig_area>0 and buffer_area==0:
            if geom.geom_type=='Polygon':
                new_geom=stitch_geom
                corr_method='"Stiching Points."'
            else:
                return non_viable

    final_area = area_me(new_geom)
    final_is_valid = new_geom.is_valid
    if final_is_valid==True:
        "Adjusted Geometry is Valid."
    final_valid_ratio = round(orig_area/final_area,4)

    return [orig_area, orig_area_type, is_valid_, is_valid_reason, new_geom, final_area, final_is_valid, final_valid_ratio]


closest_distance_rdme={
    'type':float,
    'units':'km',
    'docstring':'Takes in two geometric shapely objects and calculates the closest distance (km) between them using the haversine formula.'
}
def closest_distance(point_poly_1,point_poly_2):
    
    cen1, cen2 = nearest_points(point_poly_1, point_poly_2)
    
    lon1, lat1 = list(cen1.coords)[0]
    lon2, lat2 = list(cen2.coords)[0]

    #Fixed variables:
    radearth=6371 #In km
    
    distlat=math.radians(lat2-lat1)
    distlon=math.radians(lon2-lon1)

    a=math.sin(distlat/2)*math.sin(distlat/2)+math.cos(math.radians(lat1))*math.cos(math.radians(lat2))*math.sin(distlon/2)*math.sin(distlon/2)
    c=2*math.atan2(math.sqrt(a),math.sqrt(1-a))
    d=radearth*c
    d=round(d,4)
    
    return d


# Import service request data.
sr=pd.read_csv('20220703_231225_service_requests.csv')
sr


# Convert lat lon coordinates to Point object for intersections.
sr['geometry']=sr.apply(lambda x: Point(x['longitude'], x['latitude']), axis=1)

C:\Users\ryanp\anaconda3\lib\site-packages\pandas\core\dtypes\cast.py:122: ShapelyDeprecationWarning:

The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead.


# Create a geodataframe.
sr=gpd.GeoDataFrame(sr, geometry='geometry', crs='EPSG:4326')
sr.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 941634 entries, 0 to 941633
Data columns (total 18 columns):
 #   Column                Non-Null Count   Dtype   
---  ------                --------------   -----   
 0   Unnamed: 0.1          941634 non-null  int64   
 1   Unnamed: 0            941634 non-null  int64   
 2   notification_number   941634 non-null  int64   
 3   reference_number      592920 non-null  float64 
 4   creation_timestamp    941634 non-null  object  
 5   completion_timestamp  929442 non-null  object  
 6   directorate           932199 non-null  object  
 7   department            932180 non-null  object  
 8   branch                913233 non-null  object  
 9   section               848509 non-null  object  
 10  code_group            941634 non-null  object  
 11  code                  941634 non-null  object  
 12  cause_code_group      131117 non-null  object  
 13  cause_code            131117 non-null  object  
 14  official_suburb       729221 non-null  object  
 15  latitude              729270 non-null  float64 
 16  longitude             729270 non-null  float64 
 17  geometry              941634 non-null  geometry
dtypes: float64(3), geometry(1), int64(3), object(11)
memory usage: 129.3+ MB


# Check the coordinate reference system of the geodataframe.
sr.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich


# Check for any missing information.
msno.bar(sr)

<AxesSubplot:>


# Read in hex, resolution 8 spatial file.
hex8=gpd.read_file('20220703_231253_cpt8hex_geometries.geojson', driver='GeoJSON')
hex8


# Check that all names and geometries are unique.
hex8['index'].nunique(), hex8['geometry'].nunique()

(3832, 3832)


# Run function to check and possible correct invalid geometries.
hex8[['Orig_Area', 'Orignal Geometry Type', 'Orignal Geometry Validity', 'Orignal Geometry Validity Reason', 'Final Geometry', 'Final Area', 'Final Geometry Validity', 'Final Valid Ratio']]=hex8.apply(lambda x: polyFix_it(x['geometry'], x.name), axis=1, result_type='expand')


# Check that all geometries are valid.
hex8['Orignal Geometry Validity'].value_counts()

True    3832
Name: Orignal Geometry Validity, dtype: int64


hex8.drop(columns=['Orig_Area', 'Orignal Geometry Type', 'Orignal Geometry Validity', 'Orignal Geometry Validity Reason', 'Final Geometry', 'Final Area', 'Final Geometry Validity', 'Final Valid Ratio'], inplace=True)


# Check for any missing information.
msno.bar(hex8)

<AxesSubplot:>


# Check the coordinate reference system of the geodataframe.
hex8.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich


# Check general info of the geodataframe.
hex8.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 3832 entries, 0 to 3831
Data columns (total 4 columns):
 #   Column        Non-Null Count  Dtype   
---  ------        --------------  -----   
 0   index         3832 non-null   object  
 1   centroid_lat  3832 non-null   float64 
 2   centroid_lon  3832 non-null   float64 
 3   geometry      3832 non-null   geometry
dtypes: float64(2), geometry(1), object(1)
memory usage: 119.9+ KB


hex8['Hex area (sqkm)']=hex8['geometry'].apply(lambda x: area_me(x))


hex8['Hex area (sqkm)'].describe()

count    3832.000000
mean        0.699553
std         0.002831
min         0.692531
25%         0.697670
50%         0.699456
75%         0.702060
max         0.704613
Name: Hex area (sqkm), dtype: float64


# Join the service request layer to the hex layer.
start=time.time() 
sr_hex8=gpd.sjoin(sr, hex8[['index', 'geometry']], how='left', op='within')
duration=round(time.time()-start, 4)
print('Intersections took {} seconds'.format(duration))
sr_hex8

Intersections took 15.3435 seconds


# Calculate the number of service requests without valid coordinates.
no_coords=sr_hex8[(sr_hex8['latitude'].isna())  | (sr_hex8['longitude'].isna())].shape[0]

# Check the number entries (with valid coords) that did not intersect with the level 8 H3 layer.
viabe_sr_nohex=sr_hex8[~(sr_hex8['latitude'].isna()) | ~(sr_hex8['longitude'].isna())]
viabe_sr_nohex=viabe_sr_nohex[viabe_sr_nohex['index'].isna()].shape[0]

# Calculate total join error %.
join_error=(no_coords+viabe_sr_nohex)*100/len(sr_hex8)


# Specify the error threshold.
error_thres=50


# Log information.
logging.info('Summary of joins:')
logging.info('-'*50)
logging.info('A total of {} service requests did not have valid coordinates.'.format(no_coords))
logging.info('A total of {} service requests did not have coordinates within the boundaries of Cape Town.'.format(viabe_sr_nohex))

# Add if/else for ERROR message.
if error_thres>join_error:
    logging.info('The current join error percentage {}% does not exceed that of the threshold({}%).'.format(join_error, error_thres))
else:
    logging.error('WARNING: The current join error percentage of {}% is beyond the {}% threshold.'.format(join_error, error_thres))
logging.info('-'*50)


fig = px.pie(pd.DataFrame([['Successefully joined service requests',len(sr_hex8)-(no_coords-viabe_sr_nohex)], ['Service requests with incorrect coordinates',no_coords], ['Service requests beyond bonds of Cape Town',viabe_sr_nohex]], columns=['Join verdict','Counts']),
             values='Counts', 
             names='Join verdict', 
             title='Percentage of the join verdicts for service requests to the level 8 H3 layer')
fig.show(renderer='notebook') #Allows all plotly plots to be exported with the html.


# Rename columns.
sr_hex8.rename(columns={'index':'h3_level8_index'}, inplace=True)

# Fill no lat lons (aka no joins) with the string "0".
sr_hex8['h3_level8_index']=sr_hex8['h3_level8_index'].fillna("0")


# Read in pre-joined challenge file, "sr_hex.csv" aka "20220703_231236_service_requests_8hex.csv"
sr_8hex_valid=pd.read_csv('20220703_231236_service_requests_8hex.csv')
sr_8hex_valid


# Check counts of "0" for both geodataframes.
sr_hex8[sr_hex8['h3_level8_index']=='0'].shape[0], sr_8hex_valid[sr_8hex_valid['h3_level8_index']=='0'].shape[0]

(212367, 212364)


# Check that notification numbers are all unique.
sr_hex8['notification_number'].nunique()

941634


# Generate a merged comparative dataframe to run some checks.
compare_sr_hex=sr_hex8.merge(sr_8hex_valid, on='notification_number', how='left')
compare_sr_hex


# Extract the misaligment "0" assignments.
problem_hex8=list(compare_sr_hex[(compare_sr_hex['h3_level8_index_x']=='0') & (compare_sr_hex['h3_level8_index_y']!='0')]['h3_level8_index_y'].unique())
problem_hex8

['88ad36c629fffff', '88ad361b51fffff']


# Check the problem h3 level 8 hex from the geosjon file.
hex8[hex8['index'].isin(problem_hex8)]


# Read in hex, resolution 8 spatial file.
hex8_10=gpd.read_file('20220703_231318_cpt8_10hex_geometries.geojson', driver='GeoJSON')
hex8_10


# Check if any of the misaligned "0" level 8 hex indices occur in the large file.
hex8_10[hex8_10['index'].isin(problem_hex8)]


# Check how many entries have misaligned level 8 hex indices, excluding the 2 for the '0' mismatch.
compare_sr_hex['Same Hex yes/no']=compare_sr_hex.apply(lambda x: 'yes' if x['h3_level8_index_x']==x['h3_level8_index_y'] else 'no', axis=1)
compare_sr_hex['Same Hex yes/no'].value_counts()

yes    941605
no         29
Name: Same Hex yes/no, dtype: int64


# List these unique entries and visually inspect.
compare_sr_hex[(compare_sr_hex['Same Hex yes/no']=='no') & (compare_sr_hex['h3_level8_index_x']!='0')][['notification_number', 'latitude_x', 'longitude_x','latitude_y', 'longitude_y','h3_level8_index_x', 'h3_level8_index_y']]


# Get a unique list of the mismatched level 8 hex indices from the pre-joined file.
mismatch_hex8=list(compare_sr_hex[(compare_sr_hex['Same Hex yes/no']=='no') & (compare_sr_hex['h3_level8_index_x']!='0')]['h3_level8_index_y'].unique())
mismatch_hex8

['88ad360227fffff', '88ad36d5b5fffff', '88ad361353fffff', '88ad360157fffff']


# Check that all exists in the supplied hex8 geojson file.
[x for x in mismatch_hex8 if x not in hex8['index'].tolist()]

[]


hex8[hex8['index'].isin(mismatch_hex8)]


# Check with basic intersections if the supplied pre-joined assigments work.
[
hex8['geometry'][975].contains(Point(18.610848, -34.015341)), #88ad360157fffff
hex8['geometry'][1048].contains(Point(18.512912, -33.871389)), #88ad360227fffff
hex8['geometry'][1457].contains(Point(18.610848, -34.015341)), #88ad361353fffff
hex8['geometry'][3036].contains(Point(18.817866, -34.055004)), #88ad36d5b5fffff    
]

[False, False, True, False]


# View all 3 geometries.
GeometryCollection(hex8[hex8['index'].isin(['88ad36135bfffff','88ad361353fffff'])]['geometry'].tolist()+[(Point(18.610848, -34.015341))])


# View common area of two hexes.
hex8[hex8['index']=='88ad36135bfffff']['geometry'].item().intersection(hex8[hex8['index']=='88ad361353fffff']['geometry'].item())


# Check type of common area.
type(hex8[hex8['index']=='88ad36135bfffff']['geometry'].item().intersection(hex8[hex8['index']=='88ad361353fffff']['geometry'].item()))

shapely.geometry.linestring.LineString


# Calculate shortest distance between common linestring and Point.
closest_distance(hex8[hex8['index']=='88ad36135bfffff']['geometry'].item().intersection(hex8[hex8['index']=='88ad361353fffff']['geometry'].item()), Point(18.610848, -34.015341))

0.0


# Read in provided file.
sr_hex=pd.read_csv('20220703_231236_service_requests_8hex.csv')
sr_hex.head(10)


# Diagnose data sparsity.
msno.bar(sr_hex)

<AxesSubplot:>


# Reshape the data into number of requests created, per type, per H3 level 8 hex in the last 12 months.
sr_hex_types=sr_hex.groupby(['code', 'h3_level8_index'])['notification_number'].count().reset_index()
sr_hex_types


# Calculate the total numbers of level 8 Hexs each service request type occurs in.
sr_hex_types['code'].value_counts().reset_index()


# Calculate the total number of the service requests per type of the time period provided.
sr_hex_types.groupby(['code'])['notification_number'].sum().reset_index().sort_values(by=['notification_number'], ascending=False).reset_index(drop=True)


# Calculate the total number of weekly information per service request type for the period provided.


# Confirm before proceeding that all creation dates are in 2020.
sr_hex['year']=sr_hex['creation_timestamp'].apply(lambda x: x[:4])
sr_hex['year'].value_counts()

2020    941634
Name: year, dtype: int64


# Export week number for each creation date.
sr_hex['week']=sr_hex['creation_timestamp'].apply(lambda x: int(datetime.datetime.fromisoformat(x).strftime("%V")))
sr_hex['week'].value_counts()

29    24251
45    22864
41    22758
40    22670
42    22466
37    22335
50    22258
49    22219
4     21981
3     21876
10    21844
46    21713
36    21587
43    21507
5     21167
44    20972
34    20673
8     20641
6     20630
7     20618
48    20606
38    20341
9     20321
28    20252
11    20089
12    19578
2     19571
30    19451
47    19442
32    19354
35    19209
31    17866
27    17718
33    17594
39    17371
51    17194
26    16855
24    16297
23    14554
25    14473
13    14278
52    14237
22    13591
21    10518
53    10287
20    10220
19    10092
14     9605
1      9303
15     8979
16     8868
17     8520
18     7970
Name: week, dtype: int64


sr_hex.groupby(['code'])['week'].nunique().reset_index().sort_values(by=['week'], ascending=False).head(10).reset_index(drop=True)


# Filter for the chosen service request type.
service_rtype='Sewer: Blocked/Overflow'


sr_hex_sewers=sr_hex[sr_hex['code']==service_rtype].reset_index(drop=True)
sr_hex_sewers


# Check if all level 8 hex's for the service request type have 53 values, if not add zeros.
sr_hex_sewers.groupby(['h3_level8_index'])['week'].nunique().reset_index().sort_values(by=['week'], ascending=False)


# Specify function annotations.
adf_test_rdme={
    'type':str,
    'units':'binary - Stationary or Non-Stationary',
    'docstring':'Function takes in a df subset (service request and hex), reshapes, and runs the Augmented Dickey–Fuller test (ADF test) test to calculate the p-value and defines where the dataset is a stationary time series or not.'
}

# Define function.
def adf_test(sr_hex_subset)->adf_test_rdme:
    
    sorted_notifcations=sr_hex_subset.groupby(['week'])['notification_number'].count().reset_index().sort_values(by=['week'], ascending=True)['notification_number']
    
    result=adfuller(sorted_notifcations)
    pvalue=result[1]
    if pvalue<=0.05:
        return 'Stationary'
    else:
         return 'Non-stationary'


# Check which level 8 hex's for the chosen service request has stationary, non-stationary or too few points for the ADF test.

adf_results=[]

for hex in sr_hex_sewers['h3_level8_index'].unique():
    try:
        adf_results.append([hex,  adf_test(sr_hex_sewers[sr_hex_sewers['h3_level8_index']==hex])])
    except ValueError:
        adf_results.append([hex,  'Too little points/Division by Zero'])

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log

C:\Users\ryanp\anaconda3\lib\site-packages\statsmodels\regression\linear_model.py:924: RuntimeWarning:

divide by zero encountered in log


# Store ADF results in a dataframe.

adf_results_pd=pd.DataFrame(adf_results, columns=['h3_level8_index', 'ADF result'])
adf_results_pd


# Summarise adf results.
adf_results_pd['ADF result'].value_counts()

Stationary                            833
Non-stationary                        408
Too little points/Division by Zero    261
Name: ADF result, dtype: int64


# Sumamrise the all level 8 hex's with ADF results=='Too little points/Division by Zero' service request counts.
sr_hex_sewers[sr_hex_sewers['h3_level8_index'].isin(adf_results_pd[adf_results_pd['ADF result']=='Too little points/Division by Zero']['h3_level8_index'].unique())].groupby(['h3_level8_index'])['week'].nunique().reset_index().sort_values(by=['week'], ascending=False)


# Check the service request count distribution of all level 8 hex's with ADF results=='Too little points/Division by Zero'.
sr_hex_sewers[sr_hex_sewers['h3_level8_index'].isin(adf_results_pd[adf_results_pd['ADF result']=='Too little points/Division by Zero']['h3_level8_index'].unique())].groupby(['h3_level8_index'])['week'].nunique().reset_index().describe()


# Create subsection in logging file.
logging.info(' Hex regression error section:')
logging.info('-'*50)


# Specify function annotations.
ridge_reg_rdme={
    'type':dict,
    'units':"""multiple outputs include:
    model - Regressor Model object,
    alpha - float,
    last_week - int,
    n4weeks_predict - array,
    sorted_weeks_all - dataframe,
    fig - plotly graph object
    """,
    'docstring':"""Function runs a Ridge Regressor on the level 8 hex subset and outputs a dictionary of key objects, namely:
    model - Regressor Model,
    alpha - best alpha value chosen from tuning hyperparameters,
    last_week - last week for the service request hex subset entered,
    n4weeks_predict - array of the predicted next 4 weeks,
    sorted_weeks_all - dataframe that contains the joined dataset of the input X, Y data and the predicted X, Y data ,
    fig - plotly graph showing the "sorted_weeks_all" dataset.
    """
}

# Define function.
def ridge_reg(sr_hex_subset)->ridge_reg_rdme:
    
    # Get details of subset index.
    sr=list(sr_hex_subset['code'].unique())[0]
    hex_=list(sr_hex_subset['h3_level8_index'].unique())[0]
    
    # Setup dataset.
    sorted_weeks_nots=sr_hex_subset.groupby(['week'])['notification_number'].count().reset_index().sort_values(by=['week'], ascending=True)
    data = sorted_weeks_nots.values
    X, y = data[:, :-1], data[:, -1]
    
    # Define model evaluation method.
    cv = RepeatedKFold(n_splits=3, n_repeats=3, random_state=1)
    
    # Define model.
    model = RidgeCV(alphas=arange(0, 1, 0.01), cv=cv, scoring='neg_mean_absolute_error')
    
    # Fit model.
    try:
        model.fit(X, y)
        alpha=model.alpha_ #Best fit alpha value chosen.

        # Run prediction on next 4 weeks.
        #----------------#
        # Find last week in the subset of data.
        last_week=sorted_weeks_nots['week'].max() #Find the last reported week in the data subset.

        # Make a prediction from the last week going forward.
        n4weeks=arange(last_week, last_week+4, 1)
        n4weeks_predict=model.predict(n4weeks.reshape(-1, 1))

        # Setup figure to view output of prediction.
        predict=pd.DataFrame(list(zip(n4weeks, n4weeks_predict)), columns=['week', 'notification_number'])
        predict['Data source']='Prediction'
        sorted_weeks_nots['Data source']='Actual'
        sorted_weeks_all=sorted_weeks_nots.append(predict).reset_index(drop=True)

        fig=px.line(sorted_weeks_all, x="week", y="notification_number", color='Data source')
        #----------------#
        
        # Create dictionary for export:
        ridge_dict={
            'model':model,
            'alpha':alpha,
            'last_week_of_data':last_week,
            'predicted_next_4_weeks':n4weeks_predict,
            'Merged_dataframe':sorted_weeks_all,
            'plot':fig,
                    }
        
        return ridge_dict
    except ValueError:
        
        # Export log message.
        logging.error('Level 8 hexgaon {} for service request type {}, had too few weekly samples to proceed with the ridge regression model.'.format(hex_, sr))
        
        # Create dictionary for export:
        ridge_dict={
            'model':0,
            'alpha':0,
            'last_week_of_data':0,
            'predicted_next_4_weeks':0,
            'Merged_dataframe':0,
            'plot':0,
                    }
        
        return ridge_dict


# Specify function annotations.
auto_arima_rdme={
    'type':dict,
    'units':"""multiple outputs include:
    'model' - Auto_arima Model object,
    'last_week_of_data' - int,
    'predicted_next_4_weeks' - array,
    'Merged_dataframe' - dataframe,
    'plot' - plotly graph object
    """,
    'docstring':"""Function runs a Auto ARIMA Regressor on the level 8 hex subset and outputs a dictionary of key objects, namely:
    'model' - Auto_arima Model,
    'last_week_of_data' - last week for the service request hex subset entered,
    'predicted_next_4_weeks' - array of the predicted next 4 weeks,
    'Merged_dataframe' - dataframe that contains the joined dataset of the input X, Y data and the predicted X, Y data ,
    'plot' - plotly graph showing the "sorted_weeks_all" dataset.
    """
}

# Define function.
def auto_arima(sr_hex_subset)->auto_arima_rdme:
    
    # Get details of subset index.
    sr=list(sr_hex_subset['code'].unique())[0]
    hex_=list(sr_hex_subset['h3_level8_index'].unique())[0]
    
    # Setup dataset.
    sorted_weeks_nots=sr_hex_subset.groupby(['week'])['notification_number'].count().reset_index().sort_values(by=['week'], ascending=True)
    weeks=list(sorted_weeks_nots['week'].unique())
    last_week=max(weeks)
    
    # ARIMA assumes not breaks in data, add zero values to missing weeks.
    missing_weeks=pd.DataFrame([(x,0) for x in list(np.arange(0, 53, 1)) if x not in weeks if x<=max(weeks)], columns=['week', 'notification_number'])
    
    # Append missing set to subset.
    sorted_weeks_nots_add=sorted_weeks_nots.append(missing_weeks).sort_values(by=['week'], ascending=True).reset_index(drop=True)
    
    data = sorted_weeks_nots_add.values
    X, y = data[:, :-1], data[:, -1]
    
    # Setup Model.
    model = pm.auto_arima(y, start_p=1, start_q=1,
                      test='adf',       
                      max_p=3, max_q=3, 
                      m=1,              
                      d=None,           
                      seasonal=True,  
                      start_P=0, 
                      D=0, 
                      trace=True,
                      error_action='ignore',  
                      suppress_warnings=True, 
                      stepwise=True)
    
    # Apply prediction.
    predict_weeks = 4
    prediction, conf_int = model.predict(n_periods=predict_weeks, return_conf_int=True)
    index_of_fc = np.arange(len(y), len(y)+predict_weeks)
    
    # Setup figure to view output of prediction.
    predict=pd.DataFrame(list(zip(index_of_fc, prediction)), columns=['week', 'notification_number'])
    predict['Data source']='Prediction'
    sorted_weeks_nots_add['Data source']='Actual'
    sorted_weeks_all=sorted_weeks_nots_add.append(predict).reset_index(drop=True)

    fig=px.line(sorted_weeks_all, x="week", y="notification_number", color='Data source')
    #----------------#

    # Create dictionary for export:
    auto_arima_dict={
        'model':model,
        'last_week_of_data':last_week,
        'predicted_next_4_weeks':list(prediction),
        'Merged_dataframe':sorted_weeks_all,
        'plot':fig,
                }

    return auto_arima_dict


# Ridge Regressiion test for "Stationary"
test_stat=sr_hex_sewers[sr_hex_sewers['h3_level8_index']=='88ad36d55bfffff']


ridge_reg(test_stat)['plot']

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.


# Ridge Regressiion test for "Too little points/Division by Zero": SAMPLES greater than KFOLD.
test_low_points_grtr_kfold=sr_hex_sewers[sr_hex_sewers['h3_level8_index']=='88ad3685e7fffff']


ridge_reg(test_low_points_grtr_kfold)['plot']

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.


# Ridge Regressiion test for "Too little points/Division by Zero": SAMPLES less than KFOLD.
test_low_points=sr_hex_sewers[sr_hex_sewers['h3_level8_index']=='88ad361213fffff']


ridge_reg(test_low_points)['plot']

0


# SARIMA test of non-stationary dataset.
test_non_stat=sr_hex_sewers[sr_hex_sewers['h3_level8_index']=='88ad3612a1fffff']


auto_arima(test_non_stat)['plot']

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:36: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

Performing stepwise search to minimize aic
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.07 sec
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=217.013, Time=0.00 sec
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=217.199, Time=0.01 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=216.289, Time=0.01 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=273.895, Time=0.00 sec
 ARIMA(0,0,2)(0,0,0)[0] intercept   : AIC=216.389, Time=0.01 sec
 ARIMA(1,0,2)(0,0,0)[0] intercept   : AIC=inf, Time=0.08 sec
 ARIMA(0,0,1)(0,0,0)[0]             : AIC=262.652, Time=0.01 sec

Best model:  ARIMA(0,0,1)(0,0,0)[0] intercept
Total fit time: 0.194 seconds

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:64: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.


# Run the regression model for all the hex's of the chosen service request type for the full year.
regress_results=[]

start=time.time()
for index, row in adf_results_pd.iterrows():
    if row['ADF result']=='Stationary' or row['ADF result']=='Too little points/Division by Zero':
        regress_results.append([
            row['h3_level8_index'],
            ridge_reg(sr_hex_sewers[sr_hex_sewers['h3_level8_index']==row['h3_level8_index']])['predicted_next_4_weeks'],
            'Ridge'])
    else:
        regress_results.append([
            row['h3_level8_index'],
            auto_arima(sr_hex_sewers[sr_hex_sewers['h3_level8_index']==row['h3_level8_index']])['predicted_next_4_weeks'],
            'Ridge'])
duration=round(time.time()-start, 4)
print('Processing took {} seconds'.format(duration))

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:36: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:64: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=263.532, Time=0.02 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=269.842, Time=0.00 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=261.751, Time=0.01 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=261.546, Time=0.01 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=267.851, Time=0.00 sec
 ARIMA(0,1,2)(0,0,0)[0] intercept   : AIC=263.541, Time=0.01 sec
 ARIMA(1,1,2)(0,0,0)[0] intercept   : AIC=265.371, Time=0.02 sec
 ARIMA(0,1,1)(0,0,0)[0]             : AIC=259.548, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=261.536, Time=0.01 sec
 ARIMA(0,1,2)(0,0,0)[0]             : AIC=261.543, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0]             : AIC=259.758, Time=0.01 sec
 ARIMA(1,1,2)(0,0,0)[0]             : AIC=263.374, Time=0.02 sec

Best model:  ARIMA(0,1,1)(0,0,0)[0]          
Total fit time: 0.131 seconds

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:36: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=222.104, Time=0.01 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=239.492, Time=0.00 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=222.672, Time=0.01 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=221.770, Time=0.01 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=237.507, Time=0.00 sec
 ARIMA(0,1,2)(0,0,0)[0] intercept   : AIC=221.532, Time=0.01 sec
 ARIMA(1,1,2)(0,0,0)[0] intercept   : AIC=223.526, Time=0.02 sec
 ARIMA(0,1,3)(0,0,0)[0] intercept   : AIC=223.515, Time=0.02 sec
 ARIMA(1,1,3)(0,0,0)[0] intercept   : AIC=225.417, Time=0.03 sec
 ARIMA(0,1,2)(0,0,0)[0]             : AIC=219.622, Time=0.01 sec
 ARIMA(0,1,1)(0,0,0)[0]             : AIC=219.807, Time=0.01 sec
 ARIMA(1,1,2)(0,0,0)[0]             : AIC=221.615, Time=0.02 sec
 ARIMA(0,1,3)(0,0,0)[0]             : AIC=221.604, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=220.174, Time=0.01 sec
 ARIMA(1,1,3)(0,0,0)[0]             : AIC=223.505, Time=0.03 sec

Best model:  ARIMA(0,1,2)(0,0,0)[0]          
Total fit time: 0.221 seconds

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:64: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:36: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.05 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=226.958, Time=0.00 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=217.337, Time=0.01 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.04 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=224.962, Time=0.00 sec
 ARIMA(2,1,0)(0,0,0)[0] intercept   : AIC=205.296, Time=0.01 sec
 ARIMA(3,1,0)(0,0,0)[0] intercept   : AIC=205.568, Time=0.02 sec
 ARIMA(2,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.06 sec
 ARIMA(3,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.10 sec
 ARIMA(2,1,0)(0,0,0)[0]             : AIC=203.353, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0]             : AIC=215.342, Time=0.01 sec
 ARIMA(3,1,0)(0,0,0)[0]             : AIC=203.690, Time=0.01 sec
 ARIMA(2,1,1)(0,0,0)[0]             : AIC=200.034, Time=0.03 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=198.651, Time=0.02 sec
 ARIMA(0,1,1)(0,0,0)[0]             : AIC=196.652, Time=0.01 sec
 ARIMA(0,1,2)(0,0,0)[0]             : AIC=198.651, Time=0.02 sec
 ARIMA(1,1,2)(0,0,0)[0]             : AIC=inf, Time=0.02 sec

Best model:  ARIMA(0,1,1)(0,0,0)[0]          
Total fit time: 0.422 seconds

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:64: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:36: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.05 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=350.661, Time=0.00 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=338.135, Time=0.01 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.03 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=348.663, Time=0.00 sec
 ARIMA(2,1,0)(0,0,0)[0] intercept   : AIC=326.044, Time=0.02 sec
 ARIMA(3,1,0)(0,0,0)[0] intercept   : AIC=327.829, Time=0.02 sec
 ARIMA(2,1,1)(0,0,0)[0] intercept   : AIC=327.738, Time=0.02 sec
 ARIMA(3,1,1)(0,0,0)[0] intercept   : AIC=inf, Time=0.11 sec
 ARIMA(2,1,0)(0,0,0)[0]             : AIC=324.239, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0]             : AIC=336.145, Time=0.01 sec
 ARIMA(3,1,0)(0,0,0)[0]             : AIC=325.987, Time=0.02 sec
 ARIMA(2,1,1)(0,0,0)[0]             : AIC=325.887, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=327.565, Time=0.01 sec
 ARIMA(3,1,1)(0,0,0)[0]             : AIC=inf, Time=0.08 sec

Best model:  ARIMA(2,1,0)(0,0,0)[0]          
Total fit time: 0.399 seconds

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:64: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\3603591465.py:58: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:36: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:64: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\4236441363.py:36: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.

Performing stepwise search to minimize aic
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=127.452, Time=0.02 sec
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=126.209, Time=0.01 sec
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=125.762, Time=0.01 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=125.455, Time=0.01 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=146.040, Time=0.01 sec
 ARIMA(0,0,2)(0,0,0)[0] intercept   : AIC=127.446, Time=0.01 sec
 ARIMA(1,0,2)(0,0,0)[0] intercept   : AIC=127.034, Time=0.04 sec
 ARIMA(0,0,1)(0,0,0)[0]             : AIC=146.640, Time=0.01 sec

Best model:  ARIMA(0,0,1)(0,0,0)[0] intercept
Total fit time: 0.115 seconds
Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=347.070, Time=0.02 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=359.598, Time=0.00 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=347.607, Time=0.01 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=345.082, Time=0.01 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=357.602, Time=0.00 sec
 ARIMA(0,1,2)(0,0,0)[0] intercept   : AIC=347.075, Time=0.02 sec
 ARIMA(1,1,2)(0,0,0)[0] intercept   : AIC=348.981, Time=0.03 sec
 ARIMA(0,1,1)(0,0,0)[0]             : AIC=343.085, Time=0.01 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=345.071, Time=0.01 sec
 ARIMA(0,1,2)(0,0,0)[0]             : AIC=345.077, Time=0.01 sec
 ARIMA(1,1,0)(0,0,0)[0]             : AIC=345.608, Time=0.01 sec
 ARIMA(1,1,2)(0,0,0)[0]             : AIC=346.982, Time=0.02 sec

Best model:  ARIMA(0,1,1)(0,0,0)[0]          
Total fit time: 0.152 seconds


# Specify date and name and export the file.
date_time=datetime.datetime.now().strftime("%Y%m%d_%H%M%S")
name=date_time+'_4_week_prediction_per_hex_sewer_blocked.csv'

pd.DataFrame(regress_results, columns=['h3_level8_index', '4 week predictions', 'Model used']).to_csv(name)


# Plot count of cause codes for the chosen service request type.
fig = px.bar(sr_hex_sewers['cause_code'].value_counts().reset_index().head(20), y='cause_code', x='index', text_auto='.2s',
            title="Top 20 ranked causes of 'Sewer: Blocked/Overflow' service requests")
fig.show()


# Reload the unmodified "sr_hex" file.
sr_hex=pd.read_csv('20220703_231236_service_requests_8hex.csv')


# Specify function annotations.
c19lockdown_level_rdme={
    'type':str,
    'units':'Lockdown level',
    'docstring':"""Function takes the date of a service request in 2020 and returns a lockdown level for that request.
    """
}

# Define function.
def c19lockdown_level(date)->c19lockdown_level_rdme:
    if datetime.datetime.fromisoformat('2020-03-26 00:00:00+02:00') <= date <= datetime.datetime.fromisoformat('2020-04-30 23:59:59+02:00'):
        return 'Level 5'
    
    if datetime.datetime.fromisoformat('2020-05-01 00:00:00+02:00') <= date <= datetime.datetime.fromisoformat('2020-05-31 23:59:59+02:00'):
        return 'Level 4'
    
    if datetime.datetime.fromisoformat('2020-06-01 00:00:00+02:00') <= date <= datetime.datetime.fromisoformat('2020-08-17 23:59:59+02:00'):
        return 'Level 3'
    
    if datetime.datetime.fromisoformat('2020-08-18 00:00:00+02:00') <= date <= datetime.datetime.fromisoformat('2020-09-20 23:59:59+02:00'):
        return 'Level 2'
    
    if datetime.datetime.fromisoformat('2020-09-21 00:00:00+02:00') <= date <= datetime.datetime.fromisoformat('2020-12-28 23:59:59+02:00'):
        return 'Level 1'
    
    if datetime.datetime.fromisoformat('2020-12-29 00:00:00+02:00') <= date <= datetime.datetime.fromisoformat('2020-12-31 23:59:59+02:00'):
        return 'Level 3'
    
    return "No lockdown"


# Specify function annotations.
home_office_hours_rdme={
    'type':str,
    'units':'home or work hours',
    'docstring':"""Function in a day and hour of the service request and outputs whether it was done during home or office hours.
    """
}

# Define function.
def home_office_hours(day, hour)->home_office_hours_rdme:
    if day=='Saturday' or day=='Sunday':
        return 'home hours'
    if 6 <= hour <= 17:
        return 'office hours'
    else:
        return 'home hours'


# Specify function annotations.
c19lockdown_level_rdme={
    'type':str,
    'units':'Lockdown level',
    'docstring':"""Function takes a data in 2020 and returns a lockdown level.
    """
}

# Define function.
def home_office_location(level, home_work):
    if level=='Level 5' or level=='Level 4' or home_work=='home hours':
        return 'At home'
    else:
        return 'At work'


# Specify function annotations.
dbscan_analysis_rdme={
    'type':dict,
    'units':'Varied',
    'docstring':"""Function takes in coordinates, distance, and minimum number of samples as well as...
    ...containment method ('cluster' or '') and containment geometry attributes and supplies...
    ...a with the original point set, a dataframe of the clusters and either a weighted average...
    ...point density or containment geometry point density (eg. service requests/sqkm for a hex).
    Function uses the dbscan algorithm.
    """
}

# Define function.
def dbscan_analysis(coords, dist, samples, containment, contain_geometry_name, contain_geometry_area)->dbscan_analysis_rdme:
    
    # Record total points before clustering.
    total_points=len(coords)
    
    # Run spatial clustering for coordinates.
    coords=np.array(coords)
    kms_per_radian = 6371.0088 #Radius of earth.
    epsilon = dist/ kms_per_radian #Distance between two centroids is the diag with safety factor.
    
    db = DBSCAN(eps=epsilon, min_samples=samples, algorithm='ball_tree', metric='haversine').fit(np.radians(coords))
    
    cluster_labels = db.labels_
    num_clusters = len(set(cluster_labels))
    clusters = pd.Series([coords[cluster_labels == n] for n in range(num_clusters)])
    
    # Create cluster dataframe, sort for number of points per cluster and convert clusters to Points object.
    cluster_pd=pd.DataFrame(clusters, columns=['Clusters'])
    cluster_pd['Number of points in Cluster']=cluster_pd['Clusters'].apply(lambda x: len(x))
    cluster_pd=cluster_pd.sort_values(by=['Number of points in Cluster'], ascending=False).reset_index(drop=True)
    cluster_pd['Clusters']=cluster_pd.apply(lambda x: [Point(y) for y in x['Clusters']] if x['Number of points in Cluster']!=0 else 0, axis=1)

    # Remove empty clusters.
    cluster_pd=cluster_pd[cluster_pd['Clusters']!=0].reset_index(drop=True) #Remove unfilled list, empty cluster.

    # Create a convex hull geometry to contain all points and calculate area of convex geometry.
    cluster_pd['Convex hull geometry']=cluster_pd['Clusters'].apply(lambda x: MultiPoint(x).convex_hull)
    cluster_pd['Convex hull area sqkm']=cluster_pd['Convex hull geometry'].apply(lambda x: area_me(x))

    # Remove all other geometries which are not Polygons (geometry.geom_type - 'Polygon')
    cluster_pd['Convex hull geometry type']=cluster_pd['Convex hull geometry'].apply(lambda x: x.geom_type )
    cluster_pd=cluster_pd[cluster_pd['Convex hull geometry type']=='Polygon'].reset_index(drop=True)
    
    # Remove all excessively small Convex hull areas (< 10 sqm). This should filter out for irregular polygons, not the best method.
    cluster_pd=cluster_pd[cluster_pd['Convex hull area sqkm']>1*10**(-5)].reset_index(drop=True)
    
    # if no clusters are found after the filter is applied, exit with current information.
    if containment=='yes' and len(cluster_pd)==0:
        return {
            'Original coordinates':coords,
            'clusters_pd':'No  clusters found',
            'weighted point density':total_points/contain_geometry_area,
        }
    elif containment!='yes' and len(cluster_pd)==0:
        return {
            'Original coordinates':coords,
            'clusters_pd':'No  clusters found',
            'weighted point density':'No  containment geometry',
        }   

    # Calculate cluster density.
    cluster_pd['Cluster density - Requests per sqkm']=cluster_pd.apply(lambda x: len(x['Clusters'])/x['Convex hull area sqkm'] if x['Convex hull area sqkm']>0 else 0, axis=1)

    # Calculate weighted point density contribution for clustered points.
    cluster_pd['Weighted value']=(cluster_pd['Number of points in Cluster']/total_points)*cluster_pd['Cluster density - Requests per sqkm']
    cluster_weight=cluster_pd['Weighted value'].sum()

    if containment=='yes':
        # Calculate remaining points (free point) count, hex density and weighed density contribution to total.
        remain_points=total_points-cluster_pd['Number of points in Cluster'].sum()
        remain_area=contain_geometry_area-cluster_pd['Convex hull area sqkm'].sum()
        remain_density=remain_points/remain_area
        remain_weight=(remain_points/total_points)*remain_density
        
        # Calculate weighted density for containment polygon.
        containment_density=cluster_weight+remain_weight
        
        # Assign containment polygon name to dataframe.
        cluster_pd['Containment Geometry']=contain_geometry_name
        
        return {
            'Original coordinates':coords,
            'clusters_pd':cluster_pd,
            'weighted point density':containment_density,
        }
    else:
        return {
            'Original coordinates':coords,
            'clusters_pd':cluster_pd,
            'weighted point density':cluster_weight,
        }


# Specify function annotations.
cluster_vs_hexavg_rdme={
    'type':'Dataframe',
    'units':'Dependant on Dataframe columns.',
    'docstring':"""Function exports the results of applying either the dbscan clustering or standard service requests per sqkm per hex.
    """
}

# Define function.
def cluster_vs_hexavg(dataset, mode)->cluster_vs_hexavg_rdme:
    
    # Copy dataframe:
    df=dataset.copy()
    
    # Obtain unique hex indices.
    hex_uni=df['h3_level8_index'].unique()
    
    if mode=='cluster':
        hex_density_weights=[]

        start=time.time()
        for hex_ in hex_uni:
            hex_filt_subset=df[df['h3_level8_index']==hex_].reset_index(drop=True)
            hex_density_weights.append([hex_, dbscan_analysis(hex_filt_subset['Point geometry tuple'].tolist(), 0.02, 5, 'yes', hex_, 0.699553)['weighted point density']])
        duration=round(time.time()-start, 4)
        print('Processing took {}seconds.'.format(duration))

        return pd.DataFrame(hex_density_weights, columns=['h3_level8_index', 'Weighted Service Request per sqkm'])
        
    else:
        hex_density_weights=[]
        
        start=time.time()
        for hex_ in hex_uni:
            hex_filt_subset=df[df['h3_level8_index']==hex_].reset_index(drop=True)
            hex_density_weights.append([hex_, len(hex_filt_subset)/0.699553])
        duration=round(time.time()-start, 4)
        print('Processing took {}seconds.'.format(duration))

        return pd.DataFrame(hex_density_weights, columns=['h3_level8_index', 'Weighted Service Request per sqkm'])


# Specify function annotations.
outlier_bounds_rdme={
    'type':list,
    'units':'Dependant on input list units',
    'docstring':"""Function takes in a list of values (floats/int) and outputs the min and max value, below and above which (respectively) outliers will be found.
    """
}

# Define function.
def outlier_bounds(oned_data)->outlier_bounds_rdme:
    q75,q25 = np.percentile(oned_data,[75,25])
    intr_qr = q75-q25
    max_ = q75+(1.5*intr_qr)
    min_ = q25-(1.5*intr_qr)
    
    return [min_, max_]

# Specify function annotations.
outlier_class_rdme={
    'type':str,
    'units':'yes or no',
    'docstring':"""Function takes a value and the min max value from "outlier_bounds" and defines whether the value is an outlier or not.
    """
}

# Define function.
def outlier_class(value, min_, max_)->outlier_class_rdme:
    if min_<=value<=max_:
        return 'no'
    else:
        return 'yes'


"""
These functions were sourced from https://gist.github.com/drewda/1299198 as the library jenkspy is not compatible with pyton 3.9.
These functions were also cross-reference with that found in the jenkspy github repository, 
"https://github.com/mthh/jenkspy/blob/master/jenkspy/core.py" and the jenks repository 
"https://github.com/perrygeo/jenks/blob/master/jenks.pyx"
"""

# Specify function annotations.
getJenksBreaks_rdme={
    'type':list,
    'units':'Dependant on input list',
    'docstring':"""Function takes a list of data values as well as the number of classes and applies natural breaks to output the intervals.
    """
}

# Define function.
def getJenksBreaks(dataList, numClass)->getJenksBreaks_rdme:
  dataList.sort()
  mat1 = []
  for i in range(0,len(dataList)+1):
    temp = []
    for j in range(0,numClass+1):
      temp.append(0)
    mat1.append(temp)
  mat2 = []
  for i in range(0,len(dataList)+1):
    temp = []
    for j in range(0,numClass+1):
      temp.append(0)
    mat2.append(temp)
  for i in range(1,numClass+1):
    mat1[1][i] = 1
    mat2[1][i] = 0
    for j in range(2,len(dataList)+1):
      mat2[j][i] = float('inf')
  v = 0.0
  for l in range(2,len(dataList)+1):
    s1 = 0.0
    s2 = 0.0
    w = 0.0
    for m in range(1,l+1):
      i3 = l - m + 1
      val = float(dataList[i3-1])
      s2 += val * val
      s1 += val
      w += 1
      v = s2 - (s1 * s1) / w
      i4 = i3 - 1
      if i4 != 0:
        for j in range(2,numClass+1):
          if mat2[l][j] >= (v + mat2[i4][j - 1]):
            mat1[l][j] = i3
            mat2[l][j] = v + mat2[i4][j - 1]
    mat1[l][1] = 1
    mat2[l][1] = v
  k = len(dataList)
  kclass = []
  for i in range(0,numClass+1):
    kclass.append(0)
  kclass[numClass] = float(dataList[len(dataList) - 1])
  countNum = numClass
  while countNum >= 2:#print "rank = " + str(mat1[k][countNum])
    id = int((mat1[k][countNum]) - 2)
    #print "val = " + str(dataList[id])
    kclass[countNum - 1] = dataList[id]
    k = int((mat1[k][countNum] - 1))
    countNum -= 1
  return kclass


# Specify function annotations.
getGVF_rdme={
    'type':float,
    'units':'unitless',
    'docstring':"""Function takes a list of data values as well as the number of classes and outputs the the Goodness of Variance Fit between 0 and 1.
    """
}

# Define function.
def getGVF(dataList, numClass)->getGVF_rdme:
  """
  The Goodness of Variance Fit (GVF) is found by taking the 
  difference between the squared deviations
  from the array mean (SDAM) and the squared deviations from the 
  class means (SDCM), and dividing by the SDAM
  """
  breaks = getJenksBreaks(dataList, numClass)
  dataList.sort()
  listMean = sum(dataList)/len(dataList)
#   print listMean
  SDAM = 0.0
  for i in range(0,len(dataList)):
    sqDev = (dataList[i] - listMean)**2
    SDAM += sqDev
  SDCM = 0.0
  for i in range(0,numClass):
    if breaks[i] == 0:
      classStart = 0
    else:
      classStart = dataList.index(breaks[i])
      classStart += 1
    classEnd = dataList.index(breaks[i+1])
    classList = dataList[classStart:classEnd+1]
    classMean = sum(classList)/len(classList)
#     print classMean
    preSDCM = 0.0
    for j in range(0,len(classList)):
      sqDev2 = (classList[j] - classMean)**2
      preSDCM += sqDev2
    SDCM += preSDCM
  return (SDAM - SDCM)/SDAM
  
# Specify function annotations.
classify_rdme={
    'type':int,
    'units':'unitless',
    'docstring':"""Function takes a value and the breaks outputted from "getJenksBreaks" and assigns the interval to which the value belongs.
    """
}

# Define function.
def classify(value, breaks)->classify_rdme:
  for i in range(1, len(breaks)):
    if value < breaks[i]:
      return i
  return len(breaks) - 1


# Convert to datetime object.
sr_hex['creation datetime']=sr_hex['creation_timestamp'].apply(lambda x: datetime.datetime.fromisoformat(x))

# Extract month from datetime object.
sr_hex['creation month']=sr_hex['creation datetime'].apply(lambda x: x.strftime("%B"))

# Extract week from datetime object.
sr_hex['creation week']=sr_hex['creation datetime'].apply(lambda x: x.strftime("%V"))

# Extract day from datetime object.
sr_hex['creation day']=sr_hex['creation datetime'].apply(lambda x: x.strftime("%A"))

# Extract hour from datetime object.
sr_hex['creation hour']=sr_hex['creation datetime'].apply(lambda x: x.hour)

# Assign lockdown level.
sr_hex['Lockdown Level']=sr_hex['creation datetime'].apply(lambda x: c19lockdown_level(x))

# Assign home/office hours.
sr_hex['Home or office hours']=sr_hex.apply(lambda x: home_office_hours(x['creation day'], x['creation hour']), axis=1)

# Assign whether the call was made at home or at the office.
sr_hex['Home to Office location']=sr_hex.apply(lambda x: home_office_location(x['Lockdown Level'], x['Home or office hours']), axis=1)

# Add a key for "code lat lon month"
sr_hex['code lat lon month key']=sr_hex['code'].astype(str)+'_'+sr_hex['latitude'].astype(str)+'_'+sr_hex['longitude'].astype(str)+'_'+sr_hex['creation month'].astype(str)

# Add a key for "lat lon"
sr_hex['lat lon key']=sr_hex['latitude'].astype(str)+'_'+sr_hex['longitude'].astype(str)


# Check the request count per lockdown level.
sr_hex['Lockdown Level'].value_counts(), sr_hex['Lockdown Level'].value_counts().sum()

(Level 1        291139
 No lockdown    246837
 Level 3        209982
 Level 2        100253
 Level 4         47463
 Level 5         45960
 Name: Lockdown Level, dtype: int64,
 941634)


# Check the request count per office/home hours.
sr_hex['Home or office hours'].value_counts()

office hours    703031
home hours      238603
Name: Home or office hours, dtype: int64


# Filter all service requests for only designated at home and ignoring non-assigned hex indices.
sr_hex_filt=sr_hex[(sr_hex['Home to Office location']=='At home') & (sr_hex['h3_level8_index']!='0')].reset_index(drop=True)
# sr_hex_filt.drop_duplicates(subset=['code lat lon month key'], inplace=True)
sr_hex_filt.drop_duplicates(subset=['lat lon key'], inplace=True)
sr_hex_filt['Point geometry tuple']=sr_hex_filt.apply(lambda x: (x['longitude'], x['latitude']), axis=1)


# Run calculation for Service requests per sqkm for hex.
hex_density_weights_pd=cluster_vs_hexavg(sr_hex_filt, '')
hex_density_weights_pd

Processing took 10.5851seconds.


# Create geodataframe from hex weights.
hex_dens_gpd=gpd.GeoDataFrame(
    hex_density_weights_pd.merge(hex8[['index', 'geometry']], left_on='h3_level8_index', right_on='index', how='left'),
    geometry='geometry',
    crs='EPSG:4326'
)


# Modify geodataframe to accomodate for layering in plotly map.
hex_dens_gpd_filt=hex_dens_gpd[~(hex_dens_gpd['geometry'].isna())].reset_index(drop=True)
hex_dens_gpd_filt['Layer on/off']='on'
hex_dens_gpd_filt=hex_dens_gpd_filt.append(pd.DataFrame([['No layer', 0, 'No layer', Point(-33.80011792288416,18.61286875931004), 'off']], columns=['h3_level8_index','Weighted Service Request per sqkm','index','geometry', 'Layer on/off'])).reset_index(drop=True)
hex_dens_gpd_filt

C:\Users\ryanp\AppData\Local\Temp\ipykernel_15156\409345266.py:4: FutureWarning:

The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.


# Generate plotly map.
fig = px.choropleth_mapbox(hex_dens_gpd_filt,
                           geojson=hex_dens_gpd_filt['geometry'],
                           hover_name=hex_dens_gpd_filt['h3_level8_index'],
                           locations=hex_dens_gpd_filt.index,
                           color="Weighted Service Request per sqkm",
                           center={"lat": -33.99743, "lon": 18.53369},
                           mapbox_style="open-street-map",
#                            opacity=0.5,
                           zoom=9,
                           width=1000,
                           height=900,
                           animation_frame='Layer on/off',
                           )

fig.update_layout(
    mapbox_style="white-bg",
    mapbox_layers=[
        {
            "below": 'traces',
            "sourcetype": "raster",
            "sourceattribution": "United States Geological Survey",
            "source": [
                "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}"
            ]
        }
      ])

fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})

fig.show()


# Calculate the min max bounds for the dataset to generate outlier classification.
bounds_min, bounds_max=outlier_bounds(hex_density_weights_pd['Weighted Service Request per sqkm'].tolist())
bounds_min, bounds_max

(-301.4424925631082, 529.0878603908496)


# Review the number of outliers.
hex_dens_gpd['Outlier yes/no']=hex_dens_gpd['Weighted Service Request per sqkm'].apply(lambda x: outlier_class(x, bounds_min, bounds_max))
hex_dens_gpd['Outlier yes/no'].value_counts()

no     1716
yes      60
Name: Outlier yes/no, dtype: int64


# The highest total requests/hex area entry is taken (which is an outlier).
hex_dens_gpd[hex_dens_gpd['Weighted Service Request per sqkm']==hex_dens_gpd['Weighted Service Request per sqkm'].max()]


# Apply the clustering model to the hex.
outlier_1=dbscan_analysis(
    sr_hex_filt[sr_hex_filt['h3_level8_index']=='88ad361009fffff']['Point geometry tuple'].tolist(),
    0.02,
    5,
    'yes',
    '88ad361009fffff',
    0.699553,
)

outlier_1['clusters_pd']


# Compare the standard total requests/hex area vs the actual clustering derived 'Weighted Service Request per sqkm'.
hex_dens_gpd[hex_dens_gpd['Weighted Service Request per sqkm']==hex_dens_gpd['Weighted Service Request per sqkm'].max()]['Weighted Service Request per sqkm'].item(), outlier_1['weighted point density']

(1498.0995006811493, 7669.533510709044)


# Calculate intervals using natural breaks.
breaks=getJenksBreaks(hex_dens_gpd['Weighted Service Request per sqkm'].tolist(), 2)
breaks

[0, 231.57644953277307, 1498.0995006811493]


# Check the classification counts from the natural breaks method.
jenksclass=[classify(x, breaks) for x in hex_dens_gpd['Weighted Service Request per sqkm'].tolist()]
pd.DataFrame(jenksclass)[0].value_counts()

1    1367
2     409
Name: 0, dtype: int64


# Calculate the goodness of variance fit (gvf).
getGVF(hex_dens_gpd['Weighted Service Request per sqkm'].tolist(), 2)

0.6707682183378691


# Setup kmeans model.
X = np.array(hex_dens_gpd['Weighted Service Request per sqkm'].tolist()).reshape(-1,1)
kmeans = KMeans(n_clusters=2, random_state=0).fit(X)


# Review kmeans assignments.
pd.DataFrame(kmeans.labels_)[0].value_counts()

1    1367
0     409
Name: 0, dtype: int64


# Review centroids for kmeans model clusters.
[list(x)[0] for x in kmeans.cluster_centers_]

[394.52716954494485, 67.95435619333]


# Assign categorical classes, sparsely or densely populated.
hex_dens_gpd['Population density class']=jenksclass
hex_dens_gpd['Population density class']=hex_dens_gpd['Population density class'].apply(lambda x: 'densely populated' if x==2 else 'sparsely populated')
hex_dens_gpd['Population density class'].value_counts()

sparsely populated    1367
densely populated      409
Name: Population density class, dtype: int64


# Check that all geometries are present.
hex_dens_gpd[hex_dens_gpd['geometry'].isna()]


# Remove the entries with no geometries.
hex_dens_gpd=hex_dens_gpd[~(hex_dens_gpd['geometry'].isna())].reset_index(drop=True)


# Split into two gdfs for mapping.
hex_dens_gpd_dense=hex_dens_gpd[hex_dens_gpd['Population density class']=='densely populated'].reset_index(drop=True)
hex_dens_gpd_sparse=hex_dens_gpd[hex_dens_gpd['Population density class']=='sparsely populated'].reset_index(drop=True)


#To create Terrain map
my_map = folium.Map(
    location=[-33.96258671782152, 18.613835787647165],
    zoom_start=10,
    tiles='Stamen Terrain',
    control_scale = True
)

basemaps = {
    'Esri Satellite': folium.TileLayer(
        tiles = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
        attr = 'Esri',
        name = 'Esri Satellite',
        overlay = True,
        control = True,
        control_scale = True
    ),    
}

basemaps['Esri Satellite'].add_to(my_map)

#Add lat lon coords of mouse.
formatter = "function(num) {return L.Util.formatNum(num, 3) + ' º ';};"
MousePosition(
    position="bottomright",
    separator=" | ",
    empty_string="NaN",
    lng_first=False,
    num_digits=20,
    prefix="Coordinates:",
    lat_formatter=formatter,
    lng_formatter=formatter,
).add_to(my_map)

#GPS marker plugin.
# plugins.LocateControl().add_to(my_map)

#Add measure tool 
plugins.MeasureControl(position='topright', primary_length_unit='meters', secondary_length_unit='miles', primary_area_unit='sqmeters', secondary_area_unit='acres').add_to(my_map)

<folium.plugins.measure_control.MeasureControl at 0x1fa389cbee0>


# Settings for geojson layers. 
# ----------------------------------#
# Dense layer.
tooltip_lp_dense = GeoJsonTooltip(
    fields = ['h3_level8_index', 'Weighted Service Request per sqkm', 'Outlier yes/no', 'Population density class'],
    aliases = ['h3_level8_index', 'Weighted Service Request per sqkm', 'Outlier yes/no', 'Population density class'],
    #localize = True,
    sticky = True,
    labels = True,
    style = """
    background-color: #F0EFEF;
    border: 0.5px solid black;
    border-radius: 0.5px,
    box-shadow: 0.5px;
    """,
    max_width = 800
)

map_colors_lp_dense = lambda x: {
    'fillColor': 'blue',
    'color': 'black',
    'weight': 2,
    'fillOpacity': 0.30
}

# Sparse layer.
tooltip_lp_sparse = GeoJsonTooltip(
    fields = ['h3_level8_index', 'Weighted Service Request per sqkm', 'Outlier yes/no', 'Population density class'],
    aliases = ['h3_level8_index', 'Weighted Service Request per sqkm', 'Outlier yes/no', 'Population density class'],
    #localize = True,
    sticky = True,
    labels = True,
    style = """
    background-color: #F0EFEF;
    border: 0.5px solid black;
    border-radius: 0.5px,
    box-shadow: 0.5px;
    """,
    max_width = 800
)

map_colors_lp_sparse = lambda x: {
    'fillColor': 'yellow',
    'color': 'black',
    'weight': 2,
    'fillOpacity': 0.30
}
# ----------------------------------#
# Add Polygons.
# GeoDataFrames can be added as geojson.
folium.GeoJson(
    hex_dens_gpd_dense,
    name = 'Densely Populated Level 8 Hex',
    style_function = map_colors_lp_dense,
    tooltip = tooltip_lp_dense
).add_to(my_map)

folium.GeoJson(
    hex_dens_gpd_sparse,
    name = 'Sparsely Populated Level 8 Hex',
    style_function = map_colors_lp_sparse,
    tooltip = tooltip_lp_sparse
).add_to(my_map)

# ----------------------------------#
# Add a layer control panel to the map. #NOTE LAYER CONTROL MUST BE ADDED AFTER ALL POINTS/POLYGONS ARE ADDED TO MAP.
my_map.add_child(folium.LayerControl())

	Unnamed: 0.1	Unnamed: 0	notification_number	reference_number	creation_timestamp	completion_timestamp	directorate	department	branch	section	code_group	code	cause_code_group	cause_code	official_suburb	latitude	longitude
0	0	0	400583534	9.109492e+09	2020-10-07 06:55:18+02:00	2020-10-08 15:36:35+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area Central	District: Blaauwberg	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Wear and tear	MONTAGUE GARDENS	-33.872839	18.522488
1	1	1	400555043	9.108995e+09	2020-07-09 16:08:13+02:00	2020-07-14 14:27:01+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	SOMERSET WEST	-34.078916	18.848940
2	2	2	400589145	9.109614e+09	2020-10-27 10:21:59+02:00	2020-10-28 17:48:15+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	STRAND	-34.102242	18.821116
3	3	3	400538915	9.108601e+09	2020-03-19 06:36:06+02:00	2021-03-29 20:34:19+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area North	District : Bellville	TD Customer complaint groups	Paint Markings Lines&Signs	Road Markings	Wear and tear	RAVENSMEAD	-33.920019	18.607209
4	4	4	400568554	NaN	2020-08-25 09:48:42+02:00	2020-08-31 08:41:13+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area South	District : Athlone	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Surfacing failure	CLAREMONT	-33.987400	18.453760
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
941629	941629	941629	1016508425	9.109974e+09	2020-12-31 23:49:38+02:00	2021-01-11 11:54:42+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation Water Distribution	WATER	Leak at Valve	NaN	NaN	WOODSTOCK	-33.931571	18.452159
941630	941630	941630	1016508432	9.109975e+09	2020-12-31 23:31:11+02:00	2021-01-04 11:46:28+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	SEWER	Sewer: Blocked/Overflow	General	Foreign Objects	FISANTEKRAAL	-33.783246	18.716554
941631	941631	941631	1016508434	9.109975e+09	2020-12-31 23:58:21+02:00	2021-01-01 00:01:08+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	WATER	Burst Pipe	NaN	NaN	NaN	NaN	NaN
941632	941632	941632	1016508442	9.109975e+09	2020-12-31 23:41:57+02:00	2021-01-05 15:59:24+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	WATER MANAGEMENT DEVICE	No Water WMD	NaN	NaN	WESBANK	-33.971099	18.659831
941633	941633	941633	1016508443	9.109975e+09	2020-12-31 23:54:09+02:00	2021-02-23 10:38:21+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	WATER	Leak at Water Meter / Stopcock	NaN	NaN	BIZWENI	-34.088078	18.862550

	index	centroid_lat	centroid_lon	geometry
0	88ad361801fffff	-33.859427	18.677843	POLYGON ((18.68119 -33.86330, 18.68357 -33.859...
1	88ad361803fffff	-33.855696	18.668766	POLYGON ((18.67211 -33.85957, 18.67450 -33.855...
2	88ad361805fffff	-33.855263	18.685959	POLYGON ((18.68931 -33.85914, 18.69169 -33.855...
3	88ad361807fffff	-33.851532	18.676881	POLYGON ((18.68023 -33.85541, 18.68261 -33.851...
4	88ad361809fffff	-33.867322	18.678806	POLYGON ((18.68215 -33.87120, 18.68454 -33.867...
...	...	...	...	...
3827	88ad369715fffff	-34.353404	18.479198	POLYGON ((18.48255 -34.35726, 18.48494 -34.353...
3828	88ad369717fffff	-34.349672	18.470112	POLYGON ((18.47346 -34.35353, 18.47585 -34.349...
3829	88ad369733fffff	-34.337717	18.477288	POLYGON ((18.48063 -34.34158, 18.48303 -34.337...
3830	88ad369739fffff	-34.349293	18.487330	POLYGON ((18.49068 -34.35315, 18.49307 -34.349...
3831	88ad36973bfffff	-34.345561	18.478243	POLYGON ((18.48159 -34.34942, 18.48398 -34.345...

	Unnamed: 0.1	Unnamed: 0	notification_number	reference_number	creation_timestamp	completion_timestamp	directorate	department	branch	section	code_group	code	cause_code_group	cause_code	official_suburb	latitude	longitude	geometry	index_right	index
0	0	0	400583534	9.109492e+09	2020-10-07 06:55:18+02:00	2020-10-08 15:36:35+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area Central	District: Blaauwberg	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Wear and tear	MONTAGUE GARDENS	-33.872839	18.522488	POINT (18.52249 -33.87284)	1047.0	88ad360225fffff
1	1	1	400555043	9.108995e+09	2020-07-09 16:08:13+02:00	2020-07-14 14:27:01+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	SOMERSET WEST	-34.078916	18.848940	POINT (18.84894 -34.07892)	3055.0	88ad36d5e1fffff
2	2	2	400589145	9.109614e+09	2020-10-27 10:21:59+02:00	2020-10-28 17:48:15+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	STRAND	-34.102242	18.821116	POINT (18.82112 -34.10224)	2946.0	88ad36d437fffff
3	3	3	400538915	9.108601e+09	2020-03-19 06:36:06+02:00	2021-03-29 20:34:19+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area North	District : Bellville	TD Customer complaint groups	Paint Markings Lines&Signs	Road Markings	Wear and tear	RAVENSMEAD	-33.920019	18.607209	POINT (18.60721 -33.92002)	1247.0	88ad361133fffff
4	4	4	400568554	NaN	2020-08-25 09:48:42+02:00	2020-08-31 08:41:13+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area South	District : Athlone	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Surfacing failure	CLAREMONT	-33.987400	18.453760	POINT (18.45376 -33.98740)	2530.0	88ad361709fffff
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
941629	941629	941629	1016508425	9.109974e+09	2020-12-31 23:49:38+02:00	2021-01-11 11:54:42+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation Water Distribution	WATER	Leak at Valve	NaN	NaN	WOODSTOCK	-33.931571	18.452159	POINT (18.45216 -33.93157)	2439.0	88ad361547fffff
941630	941630	941630	1016508432	9.109975e+09	2020-12-31 23:31:11+02:00	2021-01-04 11:46:28+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	SEWER	Sewer: Blocked/Overflow	General	Foreign Objects	FISANTEKRAAL	-33.783246	18.716554	POINT (18.71655 -33.78325)	542.0	88ad3656d7fffff
941631	941631	941631	1016508434	9.109975e+09	2020-12-31 23:58:21+02:00	2021-01-01 00:01:08+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	WATER	Burst Pipe	NaN	NaN	NaN	NaN	NaN	POINT EMPTY	NaN	NaN
941632	941632	941632	1016508442	9.109975e+09	2020-12-31 23:41:57+02:00	2021-01-05 15:59:24+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	WATER MANAGEMENT DEVICE	No Water WMD	NaN	NaN	WESBANK	-33.971099	18.659831	POINT (18.65983 -33.97110)	1502.0	88ad36c49bfffff
941633	941633	941633	1016508443	9.109975e+09	2020-12-31 23:54:09+02:00	2021-02-23 10:38:21+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	WATER	Leak at Water Meter / Stopcock	NaN	NaN	BIZWENI	-34.088078	18.862550	POINT (18.86255 -34.08808)	2974.0	88ad36d517fffff

	Unnamed: 0	notification_number	reference_number	creation_timestamp	completion_timestamp	directorate	department	branch	section	code_group	code	cause_code_group	cause_code	official_suburb	latitude	longitude	h3_level8_index
0	0	400583534	9.109492e+09	2020-10-07 06:55:18+02:00	2020-10-08 15:36:35+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area Central	District: Blaauwberg	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Wear and tear	MONTAGUE GARDENS	-33.872839	18.522488	88ad360225fffff
1	1	400555043	9.108995e+09	2020-07-09 16:08:13+02:00	2020-07-14 14:27:01+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	SOMERSET WEST	-34.078916	18.848940	88ad36d5e1fffff
2	2	400589145	9.109614e+09	2020-10-27 10:21:59+02:00	2020-10-28 17:48:15+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	STRAND	-34.102242	18.821116	88ad36d437fffff
3	3	400538915	9.108601e+09	2020-03-19 06:36:06+02:00	2021-03-29 20:34:19+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area North	District : Bellville	TD Customer complaint groups	Paint Markings Lines&Signs	Road Markings	Wear and tear	RAVENSMEAD	-33.920019	18.607209	88ad361133fffff
4	4	400568554	NaN	2020-08-25 09:48:42+02:00	2020-08-31 08:41:13+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area South	District : Athlone	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Surfacing failure	CLAREMONT	-33.987400	18.453760	88ad361709fffff
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
941629	941629	1016508425	9.109974e+09	2020-12-31 23:49:38+02:00	2021-01-11 11:54:42+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation Water Distribution	WATER	Leak at Valve	NaN	NaN	WOODSTOCK	-33.931571	18.452159	88ad361547fffff
941630	941630	1016508432	9.109975e+09	2020-12-31 23:31:11+02:00	2021-01-04 11:46:28+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	SEWER	Sewer: Blocked/Overflow	General	Foreign Objects	FISANTEKRAAL	-33.783246	18.716554	88ad3656d7fffff
941631	941631	1016508434	9.109975e+09	2020-12-31 23:58:21+02:00	2021-01-01 00:01:08+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	WATER	Burst Pipe	NaN	NaN	NaN	NaN	NaN	0
941632	941632	1016508442	9.109975e+09	2020-12-31 23:41:57+02:00	2021-01-05 15:59:24+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	WATER MANAGEMENT DEVICE	No Water WMD	NaN	NaN	WESBANK	-33.971099	18.659831	88ad36c49bfffff
941633	941633	1016508443	9.109975e+09	2020-12-31 23:54:09+02:00	2021-02-23 10:38:21+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	WATER	Leak at Water Meter / Stopcock	NaN	NaN	BIZWENI	-34.088078	18.862550	88ad36d517fffff

	Unnamed: 0.1	Unnamed: 0_x	notification_number	reference_number_x	creation_timestamp_x	completion_timestamp_x	directorate_x	department_x	branch_x	section_x	...	branch_y	section_y	code_group_y	code_y	cause_code_group_y	cause_code_y	official_suburb_y	latitude_y	longitude_y	h3_level8_index_y
0	0	0	400583534	9.109492e+09	2020-10-07 06:55:18+02:00	2020-10-08 15:36:35+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area Central	District: Blaauwberg	...	RIM Area Central	District: Blaauwberg	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Wear and tear	MONTAGUE GARDENS	-33.872839	18.522488	88ad360225fffff
1	1	1	400555043	9.108995e+09	2020-07-09 16:08:13+02:00	2020-07-14 14:27:01+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	...	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	SOMERSET WEST	-34.078916	18.848940	88ad36d5e1fffff
2	2	2	400589145	9.109614e+09	2020-10-27 10:21:59+02:00	2020-10-28 17:48:15+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area East	District : Somerset West	...	RIM Area East	District : Somerset West	TD Customer complaint groups	Manhole Cover/Gully Grid	Road (RCL)	Vandalism	STRAND	-34.102242	18.821116	88ad36d437fffff
3	3	3	400538915	9.108601e+09	2020-03-19 06:36:06+02:00	2021-03-29 20:34:19+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area North	District : Bellville	...	RIM Area North	District : Bellville	TD Customer complaint groups	Paint Markings Lines&Signs	Road Markings	Wear and tear	RAVENSMEAD	-33.920019	18.607209	88ad361133fffff
4	4	4	400568554	NaN	2020-08-25 09:48:42+02:00	2020-08-31 08:41:13+02:00	URBAN MOBILITY	Roads Infrastructure Management	RIM Area South	District : Athlone	...	RIM Area South	District : Athlone	TD Customer complaint groups	Pothole&Defect Road Foot Bic Way/Kerbs	Road (RCL)	Surfacing failure	CLAREMONT	-33.987400	18.453760	88ad361709fffff
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
941629	941629	941629	1016508425	9.109974e+09	2020-12-31 23:49:38+02:00	2021-01-11 11:54:42+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation Water Distribution	...	Reticulation	Reticulation Water Distribution	WATER	Leak at Valve	NaN	NaN	WOODSTOCK	-33.931571	18.452159	88ad361547fffff
941630	941630	941630	1016508432	9.109975e+09	2020-12-31 23:31:11+02:00	2021-01-04 11:46:28+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	...	Reticulation	NaN	SEWER	Sewer: Blocked/Overflow	General	Foreign Objects	FISANTEKRAAL	-33.783246	18.716554	88ad3656d7fffff
941631	941631	941631	1016508434	9.109975e+09	2020-12-31 23:58:21+02:00	2021-01-01 00:01:08+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	...	Reticulation	NaN	WATER	Burst Pipe	NaN	NaN	NaN	NaN	NaN	0
941632	941632	941632	1016508442	9.109975e+09	2020-12-31 23:41:57+02:00	2021-01-05 15:59:24+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	...	Customer Services (Water)	Meter Management	WATER MANAGEMENT DEVICE	No Water WMD	NaN	NaN	WESBANK	-33.971099	18.659831	88ad36c49bfffff
941633	941633	941633	1016508443	9.109975e+09	2020-12-31 23:54:09+02:00	2021-02-23 10:38:21+02:00	WATER AND SANITATION	Commercial Services	Customer Services (Water)	Meter Management	...	Customer Services (Water)	Meter Management	WATER	Leak at Water Meter / Stopcock	NaN	NaN	BIZWENI	-34.088078	18.862550	88ad36d517fffff

INTRODUCTION

Import dependancies

Defined upfront functions

SECTION 1: Initial Data Transformation

Join full service requests to level 8 H3 layer

SECTION 2: Predictive Analytic Tasks

SECTION 2.1: Introspection challenge

LIMITED TIME - IMPROVEMENT NOTES

SECTION 2.2: Classification challenge

	notification_number	latitude_x	longitude_x	latitude_y	longitude_y	h3_level8_index_x	h3_level8_index_y
248172	1015706215	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
248443	1015706515	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
260685	1015720316	-34.055004	18.817866	-34.055004	18.817866	88ad36d5b1fffff	88ad36d5b5fffff
271704	1015732806	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
296215	1015760130	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
334038	1015802932	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
347502	1015818134	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
348394	1015819134	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff
360592	1015833016	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
362621	1015835341	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
363158	1015835955	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
363465	1015836310	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
366451	1015839765	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
367314	1015840714	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
391924	1015868773	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
393728	1015870909	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
398552	1015876501	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
506268	1016000828	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff
506284	1016000846	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff
506353	1016000926	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff
506392	1016000973	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff
510959	1016006213	-33.871389	18.512912	-33.871389	18.512912	88ad360221fffff	88ad360227fffff
738513	1016271817	-33.819384	18.535454	-33.819384	18.535454	88ad360151fffff	88ad360157fffff
829911	1016377419	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff
830464	1016378097	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff
830797	1016378485	-34.015341	18.610848	-34.015341	18.610848	88ad36135bfffff	88ad361353fffff

	index	centroid_lat	centroid_lon	geometry	Hex area (sqkm)
975	88ad360157fffff	-33.815443	18.534616	POLYGON ((18.53796 -33.81932, 18.54034 -33.815...	0.700345
1048	88ad360227fffff	-33.867380	18.515058	POLYGON ((18.51840 -33.87126, 18.52078 -33.867...	0.699523
1457	88ad361353fffff	-34.011408	18.610119	POLYGON ((18.61347 -34.01528, 18.61585 -34.011...	0.698107
3036	88ad36d5b5fffff	-34.053525	18.822451	POLYGON ((18.82581 -34.05739, 18.82819 -34.053...	0.698744

	code	h3_level8_index	notification_number
0	AWS: Enquiry	0	5
1	AWS: Enquiry	88ad3612c3fffff	2
2	AWS: Enquiry	88ad361887fffff	1
3	AWS: Enquiry	88ad361a31fffff	1
4	AWS: Enquiry	88ad368d03fffff	1
...	...	...	...
102707	Wires Down	88ad36d5cbfffff	5
102708	Wires Down	88ad36d5cdfffff	2
102709	Wires Down	88ad36d703fffff	1
102710	Wires Down	88ad36d711fffff	2
102711	Wires Down	88ad36d717fffff	4

	index	code
0	Sewer: Blocked/Overflow	1502
1	No Power	1400
2	No Water WMD	1327
3	No Water Supply	1312
4	Broken Leading	1301
...	...	...
487	Vendor Software (Detail in Sort Field)	1
488	Assaults	1
489	Conflicting Stages of Traffic Signals	1
490	Government Accounts	1
491	Drop-Off : Operating Times	1

	code	notification_number
0	Sewer: Blocked/Overflow	134429
1	No Power	106601
2	No Water WMD	53923
3	Street Lights - All Lights Out	39359
4	Street Lights - Single Light Out	38206
...	...	...
487	Pre-Paid Vendor Reconciliation	1
488	R003 - Difference in Unit Price	1
489	Waste Service Provider Query	1
490	Drop-Off : Operating Times	1
491	Expected Payment Date Confirmation	1

	code	week
0	Wires Down	53
1	Internal : Billing / Adjustment Queries	53
2	Sewer: Pipe Broken	53
3	Sewer: Manhole Cover: Stolen/Missing	53
4	DO NOT USE: Meter Aqualoc Fault	53
5	Damaged Bin - 240L	53
6	Sewer: Manhole Cover: Damaged	53
7	Sewer: Foul Smell	53
8	Sewer: Blocked/Overflow	53
9	Street People	53

	Unnamed: 0	notification_number	reference_number	creation_timestamp	completion_timestamp	directorate	department	branch	section	code_group	code	cause_code_group	cause_code	official_suburb	latitude	longitude	h3_level8_index	year	week
0	481	1016495961	9.109952e+09	2020-12-26 20:17:23+02:00	2020-12-27 14:48:03+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	NaN	NaN	BONTEHEUWEL	-33.949790	18.547810	88ad361033fffff	2020	52
1	603	1015764939	9.108671e+09	2020-04-11 09:54:09+02:00	2020-06-26 12:50:51+02:00	HUMAN SETTLEMENTS	Home Ownership Transfer, Tenancy Management an...	Upgrades and Maintenance	Area South	SEWER	Sewer: Blocked/Overflow	NaN	NaN	LAVENDER HILL	-34.069183	18.483887	88ad368c21fffff	2020	15
2	614	1015936877	9.108978e+09	2020-07-07 05:48:39+02:00	2020-07-07 08:26:05+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	NaN	NaN	GUGULETU	-33.975780	18.572910	88ad361001fffff	2020	28
3	631	1016411914	9.109811e+09	2020-11-29 14:23:49+02:00	2020-11-29 16:34:11+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	NaN	NaN	TAFELSIG	-34.054190	18.641240	88ad368981fffff	2020	48
4	1209	1015972707	NaN	2020-07-16 14:26:40+02:00	2020-07-19 17:17:18+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	NaN	NaN	NYANGA	-33.984780	18.576290	88ad361047fffff	2020	29
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
134424	941621	1016508414	9.109974e+09	2020-12-31 22:27:13+02:00	2021-01-11 11:37:05+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	NaN	NaN	GRACELAND	-34.041223	18.674784	88ad36c6d3fffff	2020	53
134425	941624	1016508417	9.109974e+09	2020-12-31 22:55:25+02:00	2021-01-06 08:05:40+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	NaN	NaN	AVON	-33.920677	18.566713	88ad361181fffff	2020	53
134426	941626	1016508422	9.109975e+09	2020-12-31 23:10:01+02:00	2021-01-02 08:32:18+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	General	Rags	LANGA	-33.947924	18.532356	88ad3610e1fffff	2020	53
134427	941628	1016508424	9.109975e+09	2020-12-31 23:24:26+02:00	2021-01-01 08:56:05+02:00	WATER AND SANITATION	Distribution Services	Reticulation	Reticulation WW Conveyance	SEWER	Sewer: Blocked/Overflow	NaN	NaN	NaN	NaN	NaN	0	2020	53
134428	941630	1016508432	9.109975e+09	2020-12-31 23:31:11+02:00	2021-01-04 11:46:28+02:00	WATER AND SANITATION	Distribution Services	Reticulation	NaN	SEWER	Sewer: Blocked/Overflow	General	Foreign Objects	FISANTEKRAAL	-33.783246	18.716554	88ad3656d7fffff	2020	53

	h3_level8_index	week
0	0	53
1440	88ad36d55bfffff	53
393	88ad361231fffff	53
190	88ad361003fffff	53
1400	88ad36d423fffff	53
...	...	...
318	88ad361159fffff	1
306	88ad361135fffff	1
1217	88ad368ee1fffff	1
1219	88ad368eebfffff	1
1501	88ad36d73dfffff	1

	h3_level8_index	ADF result
0	88ad361033fffff	Stationary
1	88ad368c21fffff	Stationary
2	88ad361001fffff	Stationary
3	88ad368981fffff	Stationary
4	88ad361047fffff	Stationary
...	...	...
1497	88ad361469fffff	Too little points/Division by Zero
1498	88ad360f51fffff	Too little points/Division by Zero
1499	88ad368727fffff	Too little points/Division by Zero
1500	88ad361627fffff	Too little points/Division by Zero
1501	88ad36d71bfffff	Too little points/Division by Zero

	h3_level8_index	week
82	88ad361561fffff	3
245	88ad36d531fffff	3
34	88ad360a85fffff	3
168	88ad3685e7fffff	3
240	88ad36d093fffff	3
...	...	...
113	88ad361999fffff	1
114	88ad3619e5fffff	1
116	88ad361a05fffff	1
120	88ad361a61fffff	1
260	88ad36d73dfffff	1

	h3_level8_index	Weighted Service Request per sqkm
0	88ad36d43dfffff	244.441808
1	88ad361ac7fffff	480.306710
2	88ad361149fffff	314.486536
3	88ad3612adfffff	191.550890
4	88ad361195fffff	408.832497
...	...	...
1771	88ad360aebfffff	1.429484
1772	88ad361987fffff	1.429484
1773	88ad360997fffff	1.429484
1774	88ad361469fffff	1.429484
1775	88ad360a43fffff	1.429484

	h3_level8_index	Weighted Service Request per sqkm	index	geometry	Layer on/off
0	88ad36d43dfffff	244.441808	88ad36d43dfffff	POLYGON ((18.83357 -34.12038, 18.83596 -34.116...	on
1	88ad361ac7fffff	480.306710	88ad361ac7fffff	POLYGON ((18.63442 -33.90405, 18.63680 -33.900...	on
2	88ad361149fffff	314.486536	88ad361149fffff	POLYGON ((18.64306 -33.97503, 18.64545 -33.971...	on
3	88ad3612adfffff	191.550890	88ad3612adfffff	POLYGON ((18.54180 -33.99329, 18.54419 -33.989...	on
4	88ad361195fffff	408.832497	88ad361195fffff	POLYGON ((18.54943 -33.91404, 18.55181 -33.910...	on
...	...	...	...	...	...
1771	88ad361987fffff	1.429484	88ad361987fffff	POLYGON ((18.67349 -33.80012, 18.67588 -33.796...	on
1772	88ad360997fffff	1.429484	88ad360997fffff	POLYGON ((18.57035 -33.65995, 18.57274 -33.655...	on
1773	88ad361469fffff	1.429484	88ad361469fffff	POLYGON ((18.41769 -33.96456, 18.42008 -33.960...	on
1774	88ad360a43fffff	1.429484	88ad360a43fffff	POLYGON ((18.61287 -33.79768, 18.61525 -33.793...	on
1775	No layer	0.000000	No layer	POINT (-33.80012 18.61287)	off

	Clusters	Number of points in Cluster	Convex hull geometry	Convex hull area sqkm	Convex hull geometry type	Cluster density - Requests per sqkm	Weighted value	Containment Geometry
0	[POINT (18.571461800147 -33.982997044688), POI...	273	POLYGON ((18.570675633336 -33.983808696103, 18...	0.050884	Polygon	5365.104006	1397.589116	88ad361009fffff
1	[POINT (18.56984018124404 -33.98164995859197),...	26	POLYGON ((18.570390957633 -33.982097479416, 18...	0.004030	Polygon	6451.748045	160.062452	88ad361009fffff
2	[POINT (18.569209644948 -33.983722133186), POI...	17	POLYGON ((18.56931178 -33.98375842, 18.5692096...	0.001631	Polygon	10420.872559	169.040872	88ad361009fffff
3	[POINT (18.569980149838 -33.98406503121), POIN...	15	POLYGON ((18.569988852934 -33.984294933001, 18...	0.000618	Polygon	24253.833838	347.144568	88ad361009fffff
4	[POINT (18.569417484016 -33.981279684756), POI...	13	POLYGON ((18.569415430889 -33.981686873623, 18...	0.001350	Polygon	9630.872586	119.466931	88ad361009fffff
5	[POINT (18.57141885 -33.98100595), POINT (18.5...	13	POLYGON ((18.571275277126 -33.981204732334, 18...	0.001309	Polygon	9931.832128	123.200208	88ad361009fffff
6	[POINT (18.569618333632 -33.98022167754), POIN...	12	POLYGON ((18.569617102385 -33.980465990911, 18...	0.001036	Polygon	11584.310043	132.644771	88ad361009fffff
7	[POINT (18.57185341 -33.97995229), POINT (18.5...	12	POLYGON ((18.57196569884 -33.979985547424, 18....	0.001131	Polygon	10610.556421	121.494921	88ad361009fffff
8	[POINT (18.56934034 -33.98088269), POINT (18.5...	11	POLYGON ((18.569320959101 -33.981035029314, 18...	0.000894	Polygon	12309.984255	129.207850	88ad361009fffff
9	[POINT (18.571858961024 -33.981776839112), POI...	11	POLYGON ((18.571857327809 -33.982102590286, 18...	0.000859	Polygon	12805.874992	134.412810	88ad361009fffff
10	[POINT (18.568529441093 -33.982905358068), POI...	10	POLYGON ((18.568430447427 -33.983149328463, 18...	0.000952	Polygon	10505.261457	100.241044	88ad361009fffff
11	[POINT (18.57097766 -33.98453924), POINT (18.5...	9	POLYGON ((18.570916725783 -33.984623333347, 18...	0.000641	Polygon	14032.496564	120.508081	88ad361009fffff
12	[POINT (18.570283567089 -33.983893527343), POI...	8	POLYGON ((18.57028755 -33.98390532, 18.5702835...	0.000103	Polygon	77928.560023	594.874504	88ad361009fffff
13	[POINT (18.569615206568 -33.983430257035), POI...	8	POLYGON ((18.569796204042 -33.983724185347, 18...	0.000557	Polygon	14363.946549	109.648447	88ad361009fffff
14	[POINT (18.57013932 -33.97909243), POINT (18.5...	7	POLYGON ((18.57002300241 -33.979261181676, 18....	0.000249	Polygon	28061.204352	187.431708	88ad361009fffff
15	[POINT (18.57066965 -33.98159906), POINT (18.5...	7	POLYGON ((18.570685869329 -33.981772751864, 18...	0.000267	Polygon	26229.218185	175.195160	88ad361009fffff
16	[POINT (18.569002917064 -33.982854323708), POI...	6	POLYGON ((18.569002917064 -33.982854323708, 18...	0.000288	Polygon	20838.587878	119.304892	88ad361009fffff
17	[POINT (18.571660177777 -33.982427660994), POI...	6	POLYGON ((18.571659769276 -33.982509098778, 18...	0.000251	Polygon	23914.597401	136.915634	88ad361009fffff
18	[POINT (18.57410007816404 -33.977359984589974)...	6	POLYGON ((18.57425815 -33.97742798, 18.5740943...	0.000052	Polygon	116308.726285	665.889654	88ad361009fffff
19	[POINT (18.572641022286 -33.981779557731), POI...	6	POLYGON ((18.572641022286 -33.981779557731, 18...	0.000271	Polygon	22153.137135	126.830938	88ad361009fffff
20	[POINT (18.571165660141 -33.983566087323), POI...	6	POLYGON ((18.57145853051 -33.983648546854, 18....	0.000245	Polygon	24447.576637	139.967042	88ad361009fffff
21	[POINT (18.570010177092 -33.980060168941), POI...	5	POLYGON ((18.569911191032 -33.980304140698, 18...	0.000164	Polygon	30559.489057	145.799089	88ad361009fffff
22	[POINT (18.567173612602 -33.980375981193), POI...	5	POLYGON ((18.567173612602 -33.980375981193, 18...	0.000153	Polygon	32671.106604	155.873600	88ad361009fffff
23	[POINT (18.57131532 -33.98009838), POINT (18.5...	5	POLYGON ((18.57151338 -33.98018169, 18.5711828...	0.000117	Polygon	42834.989699	204.365409	88ad361009fffff
24	[POINT (18.56823961 -33.97754163), POINT (18.5...	5	POLYGON ((18.568164761546 -33.977691968643, 18...	0.000100	Polygon	50022.720058	238.658016	88ad361009fffff
25	[POINT (18.568143759624 -33.981845295552), POI...	5	POLYGON ((18.568143759624 -33.981845295552, 18...	0.000095	Polygon	52863.600052	252.211832	88ad361009fffff
26	[POINT (18.57059803 -33.97982102), POINT (18.5...	5	POLYGON ((18.57067581 -33.97984937, 18.5704024...	0.000081	Polygon	61529.032347	293.554544	88ad361009fffff
27	[POINT (18.570097273657 -33.982177892917), POI...	5	POLYGON ((18.569999105489 -33.982258989107, 18...	0.000128	Polygon	39112.743788	186.606602	88ad361009fffff
28	[POINT (18.57205081 -33.9820766), POINT (18.57...	4	POLYGON ((18.571953861185 -33.982347243761, 18...	0.000098	Polygon	40849.417130	155.913806	88ad361009fffff
29	[POINT (18.56872454759 -33.982987481245), POIN...	3	POLYGON ((18.56872372497 -33.983150356735, 18....	0.000082	Polygon	36671.400917	104.975384	88ad361009fffff
30	[POINT (18.570490354069 -33.981772069569), POI...	3	POLYGON ((18.570588111698 -33.981772410756, 18...	0.000041	Polygon	73342.722999	209.950543	88ad361009fffff

	week
count	261.000000
mean	1.605364
std	0.755372
min	1.000000
25%	1.000000
50%	1.000000
75%	2.000000
max	3.000000